<!DOCTYPE html>
<HTML lang = "en">
<HEAD>
  <meta charset="UTF-8"/>
  <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
  <title>Ordinary Differential Equations, Applications and Discretizations</title>
  

  <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
      tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]},
      TeX: { equationNumbers: { autoNumber: "AMS" } }
    });
  </script>

  <script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
  </script>

  
<style>
pre.hljl {
    border: 1px solid #ccc;
    margin: 5px;
    padding: 5px;
    overflow-x: auto;
    color: rgb(68,68,68); background-color: rgb(251,251,251); }
pre.hljl > span.hljl-t { }
pre.hljl > span.hljl-w { }
pre.hljl > span.hljl-e { }
pre.hljl > span.hljl-eB { }
pre.hljl > span.hljl-o { }
pre.hljl > span.hljl-k { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kc { color: rgb(59,151,46); font-style: italic; }
pre.hljl > span.hljl-kd { color: rgb(214,102,97); font-style: italic; }
pre.hljl > span.hljl-kn { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kp { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kr { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kt { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-n { }
pre.hljl > span.hljl-na { }
pre.hljl > span.hljl-nb { }
pre.hljl > span.hljl-nbp { }
pre.hljl > span.hljl-nc { }
pre.hljl > span.hljl-ncB { }
pre.hljl > span.hljl-nd { color: rgb(214,102,97); }
pre.hljl > span.hljl-ne { }
pre.hljl > span.hljl-neB { }
pre.hljl > span.hljl-nf { color: rgb(66,102,213); }
pre.hljl > span.hljl-nfm { color: rgb(66,102,213); }
pre.hljl > span.hljl-np { }
pre.hljl > span.hljl-nl { }
pre.hljl > span.hljl-nn { }
pre.hljl > span.hljl-no { }
pre.hljl > span.hljl-nt { }
pre.hljl > span.hljl-nv { }
pre.hljl > span.hljl-nvc { }
pre.hljl > span.hljl-nvg { }
pre.hljl > span.hljl-nvi { }
pre.hljl > span.hljl-nvm { }
pre.hljl > span.hljl-l { }
pre.hljl > span.hljl-ld { color: rgb(148,91,176); font-style: italic; }
pre.hljl > span.hljl-s { color: rgb(201,61,57); }
pre.hljl > span.hljl-sa { color: rgb(201,61,57); }
pre.hljl > span.hljl-sb { color: rgb(201,61,57); }
pre.hljl > span.hljl-sc { color: rgb(201,61,57); }
pre.hljl > span.hljl-sd { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdB { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdC { color: rgb(201,61,57); }
pre.hljl > span.hljl-se { color: rgb(59,151,46); }
pre.hljl > span.hljl-sh { color: rgb(201,61,57); }
pre.hljl > span.hljl-si { }
pre.hljl > span.hljl-so { color: rgb(201,61,57); }
pre.hljl > span.hljl-sr { color: rgb(201,61,57); }
pre.hljl > span.hljl-ss { color: rgb(201,61,57); }
pre.hljl > span.hljl-ssB { color: rgb(201,61,57); }
pre.hljl > span.hljl-nB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nbB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nfB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nh { color: rgb(59,151,46); }
pre.hljl > span.hljl-ni { color: rgb(59,151,46); }
pre.hljl > span.hljl-nil { color: rgb(59,151,46); }
pre.hljl > span.hljl-noB { color: rgb(59,151,46); }
pre.hljl > span.hljl-oB { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-ow { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-p { }
pre.hljl > span.hljl-c { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-ch { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cm { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cp { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cpB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cs { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-csB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-g { }
pre.hljl > span.hljl-gd { }
pre.hljl > span.hljl-ge { }
pre.hljl > span.hljl-geB { }
pre.hljl > span.hljl-gh { }
pre.hljl > span.hljl-gi { }
pre.hljl > span.hljl-go { }
pre.hljl > span.hljl-gp { }
pre.hljl > span.hljl-gs { }
pre.hljl > span.hljl-gsB { }
pre.hljl > span.hljl-gt { }
</style>



  <style type="text/css">
  @font-face {
  font-style: normal;
  font-weight: 300;
}
@font-face {
  font-style: normal;
  font-weight: 400;
}
@font-face {
  font-style: normal;
  font-weight: 600;
}
html {
  font-family: sans-serif; /* 1 */
  -ms-text-size-adjust: 100%; /* 2 */
  -webkit-text-size-adjust: 100%; /* 2 */
}
body {
  margin: 0;
}
article,
aside,
details,
figcaption,
figure,
footer,
header,
hgroup,
main,
menu,
nav,
section,
summary {
  display: block;
}
audio,
canvas,
progress,
video {
  display: inline-block; /* 1 */
  vertical-align: baseline; /* 2 */
}
audio:not([controls]) {
  display: none;
  height: 0;
}
[hidden],
template {
  display: none;
}
a:active,
a:hover {
  outline: 0;
}
abbr[title] {
  border-bottom: 1px dotted;
}
b,
strong {
  font-weight: bold;
}
dfn {
  font-style: italic;
}
h1 {
  font-size: 2em;
  margin: 0.67em 0;
}
mark {
  background: #ff0;
  color: #000;
}
small {
  font-size: 80%;
}
sub,
sup {
  font-size: 75%;
  line-height: 0;
  position: relative;
  vertical-align: baseline;
}
sup {
  top: -0.5em;
}
sub {
  bottom: -0.25em;
}
img {
  border: 0;
}
svg:not(:root) {
  overflow: hidden;
}
figure {
  margin: 1em 40px;
}
hr {
  -moz-box-sizing: content-box;
  box-sizing: content-box;
  height: 0;
}
pre {
  overflow: auto;
}
code,
kbd,
pre,
samp {
  font-family: monospace, monospace;
  font-size: 1em;
}
button,
input,
optgroup,
select,
textarea {
  color: inherit; /* 1 */
  font: inherit; /* 2 */
  margin: 0; /* 3 */
}
button {
  overflow: visible;
}
button,
select {
  text-transform: none;
}
button,
html input[type="button"], /* 1 */
input[type="reset"],
input[type="submit"] {
  -webkit-appearance: button; /* 2 */
  cursor: pointer; /* 3 */
}
button[disabled],
html input[disabled] {
  cursor: default;
}
button::-moz-focus-inner,
input::-moz-focus-inner {
  border: 0;
  padding: 0;
}
input {
  line-height: normal;
}
input[type="checkbox"],
input[type="radio"] {
  box-sizing: border-box; /* 1 */
  padding: 0; /* 2 */
}
input[type="number"]::-webkit-inner-spin-button,
input[type="number"]::-webkit-outer-spin-button {
  height: auto;
}
input[type="search"] {
  -webkit-appearance: textfield; /* 1 */
  -moz-box-sizing: content-box;
  -webkit-box-sizing: content-box; /* 2 */
  box-sizing: content-box;
}
input[type="search"]::-webkit-search-cancel-button,
input[type="search"]::-webkit-search-decoration {
  -webkit-appearance: none;
}
fieldset {
  border: 1px solid #c0c0c0;
  margin: 0 2px;
  padding: 0.35em 0.625em 0.75em;
}
legend {
  border: 0; /* 1 */
  padding: 0; /* 2 */
}
textarea {
  overflow: auto;
}
optgroup {
  font-weight: bold;
}
table {
  font-family: monospace, monospace;
  font-size : 0.8em;
  border-collapse: collapse;
  border-spacing: 0;
}
td,
th {
  padding: 0;
}
thead th {
    border-bottom: 1px solid black;
    background-color: white;
}
tr:nth-child(odd){
  background-color: rgb(248,248,248);
}


/*
* Skeleton V2.0.4
* Copyright 2014, Dave Gamache
* www.getskeleton.com
* Free to use under the MIT license.
* http://www.opensource.org/licenses/mit-license.php
* 12/29/2014
*/
.container {
  position: relative;
  width: 100%;
  max-width: 960px;
  margin: 0 auto;
  padding: 0 20px;
  box-sizing: border-box; }
.column,
.columns {
  width: 100%;
  float: left;
  box-sizing: border-box; }
@media (min-width: 400px) {
  .container {
    width: 85%;
    padding: 0; }
}
@media (min-width: 550px) {
  .container {
    width: 80%; }
  .column,
  .columns {
    margin-left: 4%; }
  .column:first-child,
  .columns:first-child {
    margin-left: 0; }

  .one.column,
  .one.columns                    { width: 4.66666666667%; }
  .two.columns                    { width: 13.3333333333%; }
  .three.columns                  { width: 22%;            }
  .four.columns                   { width: 30.6666666667%; }
  .five.columns                   { width: 39.3333333333%; }
  .six.columns                    { width: 48%;            }
  .seven.columns                  { width: 56.6666666667%; }
  .eight.columns                  { width: 65.3333333333%; }
  .nine.columns                   { width: 74.0%;          }
  .ten.columns                    { width: 82.6666666667%; }
  .eleven.columns                 { width: 91.3333333333%; }
  .twelve.columns                 { width: 100%; margin-left: 0; }

  .one-third.column               { width: 30.6666666667%; }
  .two-thirds.column              { width: 65.3333333333%; }

  .one-half.column                { width: 48%; }

  /* Offsets */
  .offset-by-one.column,
  .offset-by-one.columns          { margin-left: 8.66666666667%; }
  .offset-by-two.column,
  .offset-by-two.columns          { margin-left: 17.3333333333%; }
  .offset-by-three.column,
  .offset-by-three.columns        { margin-left: 26%;            }
  .offset-by-four.column,
  .offset-by-four.columns         { margin-left: 34.6666666667%; }
  .offset-by-five.column,
  .offset-by-five.columns         { margin-left: 43.3333333333%; }
  .offset-by-six.column,
  .offset-by-six.columns          { margin-left: 52%;            }
  .offset-by-seven.column,
  .offset-by-seven.columns        { margin-left: 60.6666666667%; }
  .offset-by-eight.column,
  .offset-by-eight.columns        { margin-left: 69.3333333333%; }
  .offset-by-nine.column,
  .offset-by-nine.columns         { margin-left: 78.0%;          }
  .offset-by-ten.column,
  .offset-by-ten.columns          { margin-left: 86.6666666667%; }
  .offset-by-eleven.column,
  .offset-by-eleven.columns       { margin-left: 95.3333333333%; }

  .offset-by-one-third.column,
  .offset-by-one-third.columns    { margin-left: 34.6666666667%; }
  .offset-by-two-thirds.column,
  .offset-by-two-thirds.columns   { margin-left: 69.3333333333%; }

  .offset-by-one-half.column,
  .offset-by-one-half.columns     { margin-left: 52%; }

}
html {
  font-size: 62.5%; }
body {
  font-size: 1.5em; /* currently ems cause chrome bug misinterpreting rems on body element */
  line-height: 1.6;
  font-weight: 400;
  font-family: "Raleway", "HelveticaNeue", "Helvetica Neue", Helvetica, Arial, sans-serif;
  color: #222; }
h1, h2, h3, h4, h5, h6 {
  margin-top: 0;
  margin-bottom: 2rem;
  font-weight: 300; }
h1 { font-size: 3.6rem; line-height: 1.2;  letter-spacing: -.1rem;}
h2 { font-size: 3.4rem; line-height: 1.25; letter-spacing: -.1rem; }
h3 { font-size: 3.2rem; line-height: 1.3;  letter-spacing: -.1rem; }
h4 { font-size: 2.8rem; line-height: 1.35; letter-spacing: -.08rem; }
h5 { font-size: 2.4rem; line-height: 1.5;  letter-spacing: -.05rem; }
h6 { font-size: 1.5rem; line-height: 1.6;  letter-spacing: 0; }

p {
  margin-top: 0; }
a {
  color: #1EAEDB; }
a:hover {
  color: #0FA0CE; }
.button,
button,
input[type="submit"],
input[type="reset"],
input[type="button"] {
  display: inline-block;
  height: 38px;
  padding: 0 30px;
  color: #555;
  text-align: center;
  font-size: 11px;
  font-weight: 600;
  line-height: 38px;
  letter-spacing: .1rem;
  text-transform: uppercase;
  text-decoration: none;
  white-space: nowrap;
  background-color: transparent;
  border-radius: 4px;
  border: 1px solid #bbb;
  cursor: pointer;
  box-sizing: border-box; }
.button:hover,
button:hover,
input[type="submit"]:hover,
input[type="reset"]:hover,
input[type="button"]:hover,
.button:focus,
button:focus,
input[type="submit"]:focus,
input[type="reset"]:focus,
input[type="button"]:focus {
  color: #333;
  border-color: #888;
  outline: 0; }
.button.button-primary,
button.button-primary,
input[type="submit"].button-primary,
input[type="reset"].button-primary,
input[type="button"].button-primary {
  color: #FFF;
  background-color: #33C3F0;
  border-color: #33C3F0; }
.button.button-primary:hover,
button.button-primary:hover,
input[type="submit"].button-primary:hover,
input[type="reset"].button-primary:hover,
input[type="button"].button-primary:hover,
.button.button-primary:focus,
button.button-primary:focus,
input[type="submit"].button-primary:focus,
input[type="reset"].button-primary:focus,
input[type="button"].button-primary:focus {
  color: #FFF;
  background-color: #1EAEDB;
  border-color: #1EAEDB; }
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea,
select {
  height: 38px;
  padding: 6px 10px; /* The 6px vertically centers text on FF, ignored by Webkit */
  background-color: #fff;
  border: 1px solid #D1D1D1;
  border-radius: 4px;
  box-shadow: none;
  box-sizing: border-box; }
/* Removes awkward default styles on some inputs for iOS */
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea {
  -webkit-appearance: none;
     -moz-appearance: none;
          appearance: none; }
textarea {
  min-height: 65px;
  padding-top: 6px;
  padding-bottom: 6px; }
input[type="email"]:focus,
input[type="number"]:focus,
input[type="search"]:focus,
input[type="text"]:focus,
input[type="tel"]:focus,
input[type="url"]:focus,
input[type="password"]:focus,
textarea:focus,
select:focus {
  border: 1px solid #33C3F0;
  outline: 0; }
label,
legend {
  display: block;
  margin-bottom: .5rem;
  font-weight: 600; }
fieldset {
  padding: 0;
  border-width: 0; }
input[type="checkbox"],
input[type="radio"] {
  display: inline; }
label > .label-body {
  display: inline-block;
  margin-left: .5rem;
  font-weight: normal; }
ul {
  list-style: circle; }
ol {
  list-style: decimal; }
ul ul,
ul ol,
ol ol,
ol ul {
  margin: 1.5rem 0 1.5rem 3rem;
  font-size: 90%; }
li > p {margin : 0;}
th,
td {
  padding: 12px 15px;
  text-align: left;
  border-bottom: 1px solid #E1E1E1; }
th:first-child,
td:first-child {
  padding-left: 0; }
th:last-child,
td:last-child {
  padding-right: 0; }
button,
.button {
  margin-bottom: 1rem; }
input,
textarea,
select,
fieldset {
  margin-bottom: 1.5rem; }
pre,
blockquote,
dl,
figure,
table,
p,
ul,
ol,
form {
  margin-bottom: 1.0rem; }
.u-full-width {
  width: 100%;
  box-sizing: border-box; }
.u-max-full-width {
  max-width: 100%;
  box-sizing: border-box; }
.u-pull-right {
  float: right; }
.u-pull-left {
  float: left; }
hr {
  margin-top: 3rem;
  margin-bottom: 3.5rem;
  border-width: 0;
  border-top: 1px solid #E1E1E1; }
.container:after,
.row:after,
.u-cf {
  content: "";
  display: table;
  clear: both; }

pre {
  display: block;
  padding: 9.5px;
  margin: 0 0 10px;
  font-size: 13px;
  line-height: 1.42857143;
  word-break: break-all;
  word-wrap: break-word;
  border: 1px solid #ccc;
  border-radius: 4px;
}

pre.hljl {
  margin: 0 0 10px;
  display: block;
  background: #f5f5f5;
  border-radius: 4px;
  padding : 5px;
}

pre.output {
  background: #ffffff;
}

pre.code {
  background: #ffffff;
}

pre.julia-error {
  color : red
}

code,
kbd,
pre,
samp {
  font-family: Menlo, Monaco, Consolas, "Courier New", monospace;
  font-size: 0.9em;
}


@media (min-width: 400px) {}
@media (min-width: 550px) {}
@media (min-width: 750px) {}
@media (min-width: 1000px) {}
@media (min-width: 1200px) {}

h1.title {margin-top : 20px}
img {max-width : 100%}
div.title {text-align: center;}

  </style>
</HEAD>

<BODY>
  <div class ="container">
    <div class = "row">
      <div class = "col-md-12 twelve columns">
        <div class="title">
          <h1 class="title">Ordinary Differential Equations, Applications and Discretizations</h1>
          <h5>Chris Rackauckas</h5>
          <h5>October 1st, 2020</h5>
        </div>

        <p>Now that we have a sense of parallelism, let&#39;s return back to our thread on scientific machine learning to start constructing parallel algorithms for integration of scientific models. We previously introduced discrete dynamical systems and their asymtopic behavior. However, many physical systems are not discrete and are in fact continuous. In this discussion we will understand how to numerically compute ordinary differential equations by transforming them into discrete dynamical systems, and use this to come up with simulation techniques for physical systems.</p>
<h2>What is an Ordinary Differential Equation?</h2>
<p>An ordinary differential equation is an equation defined by a relationship on the derivative. In its general form we have that</p>
<p class="math">\[
u' = f(u,p,t)
\]</p>
<p>describes the evolution of some variable <span class="math">$u(t)$</span> which we would like to solve for. In its simplest sense, the solution to the ordinary differential equation is just the integral, since by taking the integral of both sides and applying the Fundamental Theorem of Calculus we have that</p>
<p class="math">\[
u = \int_{t_0}^{t_f} f(u,p,t)dt
\]</p>
<p>The difficulty of this equation is that the variable <span class="math">$u(t)$</span> is unknown and dependent on <span class="math">$t$</span>, meaning that the integral cannot readily be solved by simple calculus. In fact, in almost all cases there exists no analytical solution for <span class="math">$u$</span> which is readily available. However, we can understand the behavior by looking at some simple cases.</p>
<h2>Solving Ordinary Differential Equations in Julia</h2>
<p>To solve an ordinary differential equation in Julia, one can use the <a href="http://docs.juliadiffeq.org/latest/">DifferentialEquations.jl</a> package to define the differential equation you&#39;d like to solve. Let&#39;s say we want to solve the Lorenz equations:</p>
<p class="math">\[
\begin{align}
\frac{dx}{dt} &= σ(y-x) \\
\frac{dy}{dt} &= x(ρ-z) - y \\
\frac{dz}{dt} &= xy - βz \\
\end{align}
\]</p>
<p>which was the system used in our investigation of discrete dynamics. The first thing we need to do is give it this differential equation. We can either write it in an in-place form <code>f&#40;du,u,p,t&#41;</code> or an out-of-place form <code>f&#40;u,p,t&#41;</code>. Let&#39;s write it in the in-place form:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>lorenz</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span><span class='hljl-t'>
 </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>])</span><span class='hljl-t'>
 </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-oB'>-</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>])</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span>
</pre>


<pre class="output">
lorenz &#40;generic function with 1 method&#41;
</pre>


<p><strong>Question: How could I maybe speed this up a little?</strong></p>
<p>Next we give an <em>initial condition</em>. Here, this is a vector of equations, so our initial condition has to be a vector. Let&#39;s choose the following initial condition:</p>


<pre class='hljl'>
<span class='hljl-n'>u0</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>]</span>
</pre>


<pre class="output">
3-element Array&#123;Float64,1&#125;:
 1.0
 0.0
 0.0
</pre>


<p>Notice that I made sure to use <code>Float64</code> values in the initial condition. The Julia library&#39;s functions are generic and internally use the corresponding types that you give it. Integer types do not bode well for continuous problems.</p>
<p>Next, we have to tell it the timespan to solve on. Here, let&#39;s some from time 0 to 100. This means that we would use:</p>


<pre class='hljl'>
<span class='hljl-n'>tspan</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>100.0</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
&#40;0.0, 100.0&#41;
</pre>


<p>Now we need to define our parameters. We will use the same ones as from our discrete dynamical system investigation.</p>


<pre class='hljl'>
<span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>10.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>28.0</span><span class='hljl-p'>,</span><span class='hljl-ni'>8</span><span class='hljl-oB'>/</span><span class='hljl-ni'>3</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
&#40;10.0, 28.0, 2.6666666666666665&#41;
</pre>


<p>These describe an <code>ODEProblem</code>. Let&#39;s bring in DifferentialEquations.jl and define the ODE:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>DifferentialEquations</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>lorenz</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,</span><span class='hljl-n'>tspan</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
ODEProblem with uType Array&#123;Float64,1&#125; and tType Float64. In-place: true
timespan: &#40;0.0, 100.0&#41;
u0: &#91;1.0, 0.0, 0.0&#93;
</pre>


<p>Now we can solve it by calling <code>solve</code>:</p>


<pre class='hljl'>
<span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
retcode: Success
Interpolation: automatic order switching interpolation
t: 1294-element Array&#123;Float64,1&#125;:
   0.0
   3.5678604836301404e-5
   0.0003924646531993154
   0.0032624077544510573
   0.009058075635317072
   0.01695646895607931
   0.0276899566248403
   0.041856345938267966
   0.06024040228733675
   0.08368539694547242
   ⋮
  99.39403070915297
  99.47001147494375
  99.54379656909015
  99.614651558349
  99.69093823148101
  99.78733023233721
  99.86114450046736
  99.96115759510786
 100.0
u: 1294-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;1.0, 0.0, 0.0&#93;
 &#91;0.9996434557625105, 0.0009988049817849058, 1.781434788799208e-8&#93;
 &#91;0.9961045497425811, 0.010965399721242457, 2.146955365838907e-6&#93;
 &#91;0.9693591634199452, 0.08977060667778931, 0.0001438018342266937&#93;
 &#91;0.9242043615038835, 0.24228912482984957, 0.0010461623302512404&#93;
 &#91;0.8800455868998046, 0.43873645009348244, 0.0034242593451028745&#93;
 &#91;0.8483309877783048, 0.69156288756671, 0.008487623500490047&#93;
 &#91;0.8495036595681027, 1.0145425335433382, 0.01821208597613427&#93;
 &#91;0.9139069079152129, 1.4425597546855036, 0.03669381053327124&#93;
 &#91;1.0888636764765296, 2.052326153029042, 0.07402570506414284&#93;
 ⋮
 &#91;12.999157033749652, 14.10699925404482, 31.74244844521858&#93;
 &#91;11.646131422021162, 7.2855792145502845, 35.365000488215486&#93;
 &#91;7.777555445486692, 2.5166095828739574, 32.030953593541675&#93;
 &#91;4.739741627223412, 1.5919220588229062, 27.249779003951755&#93;
 &#91;3.2351668945618774, 2.3121727966182695, 22.724936101772805&#93;
 &#91;3.310411964698304, 4.28106626744641, 18.435441144016366&#93;
 &#91;4.527117863517627, 6.895878639772805, 16.58544600757436&#93;
 &#91;8.043672261487556, 12.711555298531689, 18.12537420595938&#93;
 &#91;9.97537965430362, 15.143884806010783, 21.00643286956427&#93;
</pre>


<p>To see what the solution looks like, we can call <code>plot</code>:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>Plots</span><span class='hljl-t'>
</span><span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>)</span>
</pre>


<img src=""  />

<p>We can also plot phase space diagrams by telling it which <code>vars</code> to compare on which axis. Let&#39;s plot this in the <code>&#40;x,y,z&#41;</code> plane:</p>


<pre class='hljl'>
<span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>,</span><span class='hljl-n'>vars</span><span class='hljl-oB'>=</span><span class='hljl-p'>(</span><span class='hljl-ni'>1</span><span class='hljl-p'>,</span><span class='hljl-ni'>2</span><span class='hljl-p'>,</span><span class='hljl-ni'>3</span><span class='hljl-p'>))</span>
</pre>


<img src=""  />

<p>Note that the sentinal to time is <code>0</code>, so we can also do <code>&#40;t,y,z&#41;</code> with:</p>


<pre class='hljl'>
<span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>,</span><span class='hljl-n'>vars</span><span class='hljl-oB'>=</span><span class='hljl-p'>(</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>2</span><span class='hljl-p'>,</span><span class='hljl-ni'>3</span><span class='hljl-p'>))</span>
</pre>


<img src=""  />

<p>The equation is continuous and therefore the solution is continuous. We can see this by checking how it is at any random time value:</p>


<pre class='hljl'>
<span class='hljl-nf'>sol</span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.5</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
3-element Array&#123;Float64,1&#125;:
  6.5036548742051625
 -8.508354680614291
 38.09199727173416
</pre>


<p>which gives the current evolution at that time point.</p>
<h2>Differential Equations from Scientific Contexts</h2>
<h3>N-Body Problems and Astronomy</h3>
<p>There are many different contexts in which differential equations show up. In fact, it&#39;s not a stretch to say that the laws in all fields of science are encoded in differential equations. The starting point for physics is Newton&#39;s laws of gravity, which define an N-body ordinary differential equation system by describing the force between two particles as:</p>
<p class="math">\[
F = G \frac{m_1m_2}{r^2}
\]</p>
<p>where <span class="math">$r^2$</span> is the Euclidian distance between the two particles. From here, we use the fact that</p>
<p class="math">\[
F = ma
\]</p>
<p>to receive differential equations in terms of the accelerations of each particle. The differential equation is a system, where we know the change in position is due to the current velocity:</p>
<p class="math">\[
x' = v
\]</p>
<p>and the change in velocity is the acceleration:</p>
<p class="math">\[
v' = F/m = G \frac{m_i}{r_i^2}
\]</p>
<p>where <span class="math">$i$</span> runs over the other particles. Thus we have a vector of position derivatives and a vector of velocity derivatives that evolve over time to give the evolving positions and velocity.</p>
<p>An example of this is the <a href="https://archimede.dm.uniba.it/~testset/report/plei.pdf">Pleiades problem</a>, which is an approximation to a 7-star chaotic system. It can be written as:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>OrdinaryDiffEq</span><span class='hljl-t'>

</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>pleiades</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-nd'>@inbounds</span><span class='hljl-t'> </span><span class='hljl-k'>begin</span><span class='hljl-t'>
  </span><span class='hljl-n'>x</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>view</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>7</span><span class='hljl-p'>)</span><span class='hljl-t'>   </span><span class='hljl-cs'># x</span><span class='hljl-t'>
  </span><span class='hljl-n'>y</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>view</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-ni'>8</span><span class='hljl-oB'>:</span><span class='hljl-ni'>14</span><span class='hljl-p'>)</span><span class='hljl-t'>  </span><span class='hljl-cs'># y</span><span class='hljl-t'>
  </span><span class='hljl-n'>v</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>view</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-ni'>15</span><span class='hljl-oB'>:</span><span class='hljl-ni'>21</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-cs'># x′</span><span class='hljl-t'>
  </span><span class='hljl-n'>w</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>view</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-ni'>22</span><span class='hljl-oB'>:</span><span class='hljl-ni'>28</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-cs'># y′</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>7</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>.=</span><span class='hljl-t'> </span><span class='hljl-n'>v</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>8</span><span class='hljl-oB'>:</span><span class='hljl-ni'>14</span><span class='hljl-p'>]</span><span class='hljl-oB'>.=</span><span class='hljl-t'> </span><span class='hljl-n'>w</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>15</span><span class='hljl-oB'>:</span><span class='hljl-ni'>28</span><span class='hljl-t'>
    </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>zero</span><span class='hljl-p'>(</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>])</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-oB'>=</span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>7</span><span class='hljl-p'>,</span><span class='hljl-n'>j</span><span class='hljl-oB'>=</span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>7</span><span class='hljl-t'>
    </span><span class='hljl-k'>if</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-oB'>!=</span><span class='hljl-t'> </span><span class='hljl-n'>j</span><span class='hljl-t'>
      </span><span class='hljl-n'>r</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>((</span><span class='hljl-n'>x</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>]</span><span class='hljl-oB'>-</span><span class='hljl-n'>x</span><span class='hljl-p'>[</span><span class='hljl-n'>j</span><span class='hljl-p'>])</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-n'>j</span><span class='hljl-p'>])</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-p'>)</span><span class='hljl-oB'>^</span><span class='hljl-p'>(</span><span class='hljl-ni'>3</span><span class='hljl-oB'>/</span><span class='hljl-ni'>2</span><span class='hljl-p'>)</span><span class='hljl-t'>
      </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>14</span><span class='hljl-oB'>+</span><span class='hljl-n'>i</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-n'>j</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-p'>[</span><span class='hljl-n'>j</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>x</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>])</span><span class='hljl-oB'>/</span><span class='hljl-n'>r</span><span class='hljl-t'>
      </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>21</span><span class='hljl-oB'>+</span><span class='hljl-n'>i</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-n'>j</span><span class='hljl-oB'>*</span><span class='hljl-p'>(</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-n'>j</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>])</span><span class='hljl-oB'>/</span><span class='hljl-n'>r</span><span class='hljl-t'>
    </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>tspan</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>pleiades</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>4.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>4.0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.75</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>1.5</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>1.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>],</span><span class='hljl-n'>tspan</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
ODEProblem with uType Array&#123;Float64,1&#125; and tType Float64. In-place: true
timespan: &#40;0.0, 3.0&#41;
u0: &#91;3.0, 3.0, -1.0, -3.0, 2.0, -2.0, 2.0, 3.0, -3.0, 2.0  …  0.0, 1.75, -1
.5, 0.0, 0.0, 0.0, -1.25, 1.0, 0.0, 0.0&#93;
</pre>


<p>where we assume <span class="math">$m_i = i$</span>. When we solve this equation we receive the following:</p>


<pre class='hljl'>
<span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>,</span><span class='hljl-nf'>Vern8</span><span class='hljl-p'>(),</span><span class='hljl-n'>abstol</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>1e-10</span><span class='hljl-p'>,</span><span class='hljl-n'>reltol</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>1e-10</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>)</span>
</pre>


<img src=""  />


<pre class='hljl'>
<span class='hljl-n'>tspan</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>200.0</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>pleiades</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>4.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>4.0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.75</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>1.5</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-oB'>-</span><span class='hljl-nfB'>1.25</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>,</span><span class='hljl-ni'>0</span><span class='hljl-p'>],</span><span class='hljl-n'>tspan</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>,</span><span class='hljl-nf'>Vern8</span><span class='hljl-p'>(),</span><span class='hljl-n'>abstol</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>1e-10</span><span class='hljl-p'>,</span><span class='hljl-n'>reltol</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>1e-10</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>,</span><span class='hljl-n'>vars</span><span class='hljl-oB'>=</span><span class='hljl-p'>((</span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>7</span><span class='hljl-p'>),(</span><span class='hljl-ni'>8</span><span class='hljl-oB'>:</span><span class='hljl-ni'>14</span><span class='hljl-p'>)))</span>
</pre>


<img src=""  />

<h3>Population Ecology: Lotka-Volterra</h3>
<p>Population ecology&#39;s starting point is the Lotka-Volterra equations which describes the interactions between a predator and a prey. In this case, the prey grows at an exponential rate but has a term that reduces its population by being eaten by the predator. The predator&#39;s growth is dependent on the available food &#40;the amount of prey&#41; and has a decay rate due to old age. This model is then written as follows:</p>


<pre class='hljl'>
<span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>lotka</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-ni'>4</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>u</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>

</span><span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>1.5</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>3.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>lotka</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>],(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>10.0</span><span class='hljl-p'>),</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>)</span>
</pre>


<img src=""  />

<h3>Biochemistry: Robertson Equations</h3>
<p>Biochemical equations commonly display large separation of timescales which lead to a stiffness phonomena that will be investigated later. The classic &quot;hard&quot; equations for ODE integration thus tend to come from biology &#40;not physics&#33;&#41; due to this property. One of the standard models is the Robertson model, which can be described as:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>Sundials</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>ParameterizedFunctions</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>rober</span><span class='hljl-p'>(</span><span class='hljl-n'>du</span><span class='hljl-p'>,</span><span class='hljl-n'>u</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-n'>y₁</span><span class='hljl-p'>,</span><span class='hljl-n'>y₂</span><span class='hljl-p'>,</span><span class='hljl-n'>y₃</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>u</span><span class='hljl-t'>
  </span><span class='hljl-n'>k₁</span><span class='hljl-p'>,</span><span class='hljl-n'>k₂</span><span class='hljl-p'>,</span><span class='hljl-n'>k₃</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>k₁</span><span class='hljl-oB'>*</span><span class='hljl-n'>y₁</span><span class='hljl-oB'>+</span><span class='hljl-n'>k₃</span><span class='hljl-oB'>*</span><span class='hljl-n'>y₂</span><span class='hljl-oB'>*</span><span class='hljl-n'>y₃</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'>  </span><span class='hljl-n'>k₁</span><span class='hljl-oB'>*</span><span class='hljl-n'>y₁</span><span class='hljl-oB'>-</span><span class='hljl-n'>k₂</span><span class='hljl-oB'>*</span><span class='hljl-n'>y₂</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-oB'>-</span><span class='hljl-n'>k₃</span><span class='hljl-oB'>*</span><span class='hljl-n'>y₂</span><span class='hljl-oB'>*</span><span class='hljl-n'>y₃</span><span class='hljl-t'>
  </span><span class='hljl-n'>du</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'>  </span><span class='hljl-n'>k₂</span><span class='hljl-oB'>*</span><span class='hljl-n'>y₂</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>rober</span><span class='hljl-p'>,[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>],(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1e5</span><span class='hljl-p'>),(</span><span class='hljl-nfB'>0.04</span><span class='hljl-p'>,</span><span class='hljl-nfB'>3e7</span><span class='hljl-p'>,</span><span class='hljl-nfB'>1e4</span><span class='hljl-p'>))</span><span class='hljl-t'>
</span><span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>,</span><span class='hljl-nf'>Rosenbrock23</span><span class='hljl-p'>())</span><span class='hljl-t'>
</span><span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>)</span>
</pre>


<img src=""  />


<pre class='hljl'>
<span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>xscale</span><span class='hljl-oB'>=:</span><span class='hljl-n'>log10</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>tspan</span><span class='hljl-oB'>=</span><span class='hljl-p'>(</span><span class='hljl-nfB'>1e-6</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>1e5</span><span class='hljl-p'>),</span><span class='hljl-t'> </span><span class='hljl-n'>layout</span><span class='hljl-oB'>=</span><span class='hljl-p'>(</span><span class='hljl-ni'>3</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>))</span>
</pre>


<img src=""  />

<h3>Chemical Physics: Pollution Models</h3>
<p>Chemical reactions in physical models are also described as differential equation systems. The following is a classic model of dynamics between different species of pollutants:</p>


<pre class='hljl'>
<span class='hljl-n'>k1</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.35e0</span><span class='hljl-t'>
</span><span class='hljl-n'>k2</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.266e2</span><span class='hljl-t'>
</span><span class='hljl-n'>k3</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.123e5</span><span class='hljl-t'>
</span><span class='hljl-n'>k4</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.86e-3</span><span class='hljl-t'>
</span><span class='hljl-n'>k5</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.82e-3</span><span class='hljl-t'>
</span><span class='hljl-n'>k6</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.15e5</span><span class='hljl-t'>
</span><span class='hljl-n'>k7</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.13e-3</span><span class='hljl-t'>
</span><span class='hljl-n'>k8</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.24e5</span><span class='hljl-t'>
</span><span class='hljl-n'>k9</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.165e5</span><span class='hljl-t'>
</span><span class='hljl-n'>k10</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.9e4</span><span class='hljl-t'>
</span><span class='hljl-n'>k11</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.22e-1</span><span class='hljl-t'>
</span><span class='hljl-n'>k12</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.12e5</span><span class='hljl-t'>
</span><span class='hljl-n'>k13</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.188e1</span><span class='hljl-t'>
</span><span class='hljl-n'>k14</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.163e5</span><span class='hljl-t'>
</span><span class='hljl-n'>k15</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.48e7</span><span class='hljl-t'>
</span><span class='hljl-n'>k16</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.35e-3</span><span class='hljl-t'>
</span><span class='hljl-n'>k17</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.175e-1</span><span class='hljl-t'>
</span><span class='hljl-n'>k18</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.1e9</span><span class='hljl-t'>
</span><span class='hljl-n'>k19</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.444e12</span><span class='hljl-t'>
</span><span class='hljl-n'>k20</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.124e4</span><span class='hljl-t'>
</span><span class='hljl-n'>k21</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.21e1</span><span class='hljl-t'>
</span><span class='hljl-n'>k22</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.578e1</span><span class='hljl-t'>
</span><span class='hljl-n'>k23</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.474e-1</span><span class='hljl-t'>
</span><span class='hljl-n'>k24</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.178e4</span><span class='hljl-t'>
</span><span class='hljl-n'>k25</span><span class='hljl-oB'>=</span><span class='hljl-nfB'>.312e1</span><span class='hljl-t'>
</span><span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-n'>k1</span><span class='hljl-p'>,</span><span class='hljl-n'>k2</span><span class='hljl-p'>,</span><span class='hljl-n'>k3</span><span class='hljl-p'>,</span><span class='hljl-n'>k4</span><span class='hljl-p'>,</span><span class='hljl-n'>k5</span><span class='hljl-p'>,</span><span class='hljl-n'>k6</span><span class='hljl-p'>,</span><span class='hljl-n'>k7</span><span class='hljl-p'>,</span><span class='hljl-n'>k8</span><span class='hljl-p'>,</span><span class='hljl-n'>k9</span><span class='hljl-p'>,</span><span class='hljl-n'>k10</span><span class='hljl-p'>,</span><span class='hljl-n'>k11</span><span class='hljl-p'>,</span><span class='hljl-n'>k12</span><span class='hljl-p'>,</span><span class='hljl-n'>k13</span><span class='hljl-p'>,</span><span class='hljl-n'>k14</span><span class='hljl-p'>,</span><span class='hljl-n'>k15</span><span class='hljl-p'>,</span><span class='hljl-n'>k16</span><span class='hljl-p'>,</span><span class='hljl-n'>k17</span><span class='hljl-p'>,</span><span class='hljl-n'>k18</span><span class='hljl-p'>,</span><span class='hljl-n'>k19</span><span class='hljl-p'>,</span><span class='hljl-n'>k20</span><span class='hljl-p'>,</span><span class='hljl-n'>k21</span><span class='hljl-p'>,</span><span class='hljl-n'>k22</span><span class='hljl-p'>,</span><span class='hljl-n'>k23</span><span class='hljl-p'>,</span><span class='hljl-n'>k24</span><span class='hljl-p'>,</span><span class='hljl-n'>k25</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>dy</span><span class='hljl-p'>,</span><span class='hljl-n'>y</span><span class='hljl-p'>,</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span><span class='hljl-t'>
 </span><span class='hljl-n'>k1</span><span class='hljl-p'>,</span><span class='hljl-n'>k2</span><span class='hljl-p'>,</span><span class='hljl-n'>k3</span><span class='hljl-p'>,</span><span class='hljl-n'>k4</span><span class='hljl-p'>,</span><span class='hljl-n'>k5</span><span class='hljl-p'>,</span><span class='hljl-n'>k6</span><span class='hljl-p'>,</span><span class='hljl-n'>k7</span><span class='hljl-p'>,</span><span class='hljl-n'>k8</span><span class='hljl-p'>,</span><span class='hljl-n'>k9</span><span class='hljl-p'>,</span><span class='hljl-n'>k10</span><span class='hljl-p'>,</span><span class='hljl-n'>k11</span><span class='hljl-p'>,</span><span class='hljl-n'>k12</span><span class='hljl-p'>,</span><span class='hljl-n'>k13</span><span class='hljl-p'>,</span><span class='hljl-n'>k14</span><span class='hljl-p'>,</span><span class='hljl-n'>k15</span><span class='hljl-p'>,</span><span class='hljl-n'>k16</span><span class='hljl-p'>,</span><span class='hljl-n'>k17</span><span class='hljl-p'>,</span><span class='hljl-n'>k18</span><span class='hljl-p'>,</span><span class='hljl-n'>k19</span><span class='hljl-p'>,</span><span class='hljl-n'>k20</span><span class='hljl-p'>,</span><span class='hljl-n'>k21</span><span class='hljl-p'>,</span><span class='hljl-n'>k22</span><span class='hljl-p'>,</span><span class='hljl-n'>k23</span><span class='hljl-p'>,</span><span class='hljl-n'>k24</span><span class='hljl-p'>,</span><span class='hljl-n'>k25</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-t'>
 </span><span class='hljl-n'>r1</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k1</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r2</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k2</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>4</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r3</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k3</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>5</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r4</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k4</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>7</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r5</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k5</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>7</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r6</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k6</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>7</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>6</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r7</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k7</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>9</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r8</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k8</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>9</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>6</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r9</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k9</span><span class='hljl-t'> </span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>11</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r10</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k10</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>11</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r11</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k11</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>13</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r12</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k12</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>10</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r13</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k13</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>14</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r14</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k14</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>6</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r15</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k15</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r16</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k16</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>4</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r17</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k17</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>4</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r18</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k18</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>16</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r19</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k19</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>16</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r20</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k20</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>17</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>6</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r21</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k21</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>19</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r22</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k22</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>19</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r23</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k23</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>4</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r24</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k24</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>19</span><span class='hljl-p'>]</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'>
 </span><span class='hljl-n'>r25</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>k25</span><span class='hljl-oB'>*</span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-ni'>20</span><span class='hljl-p'>]</span><span class='hljl-t'>

 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r1</span><span class='hljl-oB'>-</span><span class='hljl-n'>r10</span><span class='hljl-oB'>-</span><span class='hljl-n'>r14</span><span class='hljl-oB'>-</span><span class='hljl-n'>r23</span><span class='hljl-oB'>-</span><span class='hljl-n'>r24</span><span class='hljl-oB'>+</span><span class='hljl-t'>
          </span><span class='hljl-n'>r2</span><span class='hljl-oB'>+</span><span class='hljl-n'>r3</span><span class='hljl-oB'>+</span><span class='hljl-n'>r9</span><span class='hljl-oB'>+</span><span class='hljl-n'>r11</span><span class='hljl-oB'>+</span><span class='hljl-n'>r12</span><span class='hljl-oB'>+</span><span class='hljl-n'>r22</span><span class='hljl-oB'>+</span><span class='hljl-n'>r25</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r2</span><span class='hljl-oB'>-</span><span class='hljl-n'>r3</span><span class='hljl-oB'>-</span><span class='hljl-n'>r9</span><span class='hljl-oB'>-</span><span class='hljl-n'>r12</span><span class='hljl-oB'>+</span><span class='hljl-n'>r1</span><span class='hljl-oB'>+</span><span class='hljl-n'>r21</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>3</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r15</span><span class='hljl-oB'>+</span><span class='hljl-n'>r1</span><span class='hljl-oB'>+</span><span class='hljl-n'>r17</span><span class='hljl-oB'>+</span><span class='hljl-n'>r19</span><span class='hljl-oB'>+</span><span class='hljl-n'>r22</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>4</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r2</span><span class='hljl-oB'>-</span><span class='hljl-n'>r16</span><span class='hljl-oB'>-</span><span class='hljl-n'>r17</span><span class='hljl-oB'>-</span><span class='hljl-n'>r23</span><span class='hljl-oB'>+</span><span class='hljl-n'>r15</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>5</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r3</span><span class='hljl-oB'>+</span><span class='hljl-n'>r4</span><span class='hljl-oB'>+</span><span class='hljl-n'>r4</span><span class='hljl-oB'>+</span><span class='hljl-n'>r6</span><span class='hljl-oB'>+</span><span class='hljl-n'>r7</span><span class='hljl-oB'>+</span><span class='hljl-n'>r13</span><span class='hljl-oB'>+</span><span class='hljl-n'>r20</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>6</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r6</span><span class='hljl-oB'>-</span><span class='hljl-n'>r8</span><span class='hljl-oB'>-</span><span class='hljl-n'>r14</span><span class='hljl-oB'>-</span><span class='hljl-n'>r20</span><span class='hljl-oB'>+</span><span class='hljl-n'>r3</span><span class='hljl-oB'>+</span><span class='hljl-n'>r18</span><span class='hljl-oB'>+</span><span class='hljl-n'>r18</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>7</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r4</span><span class='hljl-oB'>-</span><span class='hljl-n'>r5</span><span class='hljl-oB'>-</span><span class='hljl-n'>r6</span><span class='hljl-oB'>+</span><span class='hljl-n'>r13</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>8</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>r4</span><span class='hljl-oB'>+</span><span class='hljl-n'>r5</span><span class='hljl-oB'>+</span><span class='hljl-n'>r6</span><span class='hljl-oB'>+</span><span class='hljl-n'>r7</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>9</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r7</span><span class='hljl-oB'>-</span><span class='hljl-n'>r8</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>10</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r12</span><span class='hljl-oB'>+</span><span class='hljl-n'>r7</span><span class='hljl-oB'>+</span><span class='hljl-n'>r9</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>11</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r9</span><span class='hljl-oB'>-</span><span class='hljl-n'>r10</span><span class='hljl-oB'>+</span><span class='hljl-n'>r8</span><span class='hljl-oB'>+</span><span class='hljl-n'>r11</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>12</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>r9</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>13</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r11</span><span class='hljl-oB'>+</span><span class='hljl-n'>r10</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>14</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r13</span><span class='hljl-oB'>+</span><span class='hljl-n'>r12</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>15</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>r14</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>16</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r18</span><span class='hljl-oB'>-</span><span class='hljl-n'>r19</span><span class='hljl-oB'>+</span><span class='hljl-n'>r16</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>17</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r20</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>18</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>r20</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>19</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r21</span><span class='hljl-oB'>-</span><span class='hljl-n'>r22</span><span class='hljl-oB'>-</span><span class='hljl-n'>r24</span><span class='hljl-oB'>+</span><span class='hljl-n'>r23</span><span class='hljl-oB'>+</span><span class='hljl-n'>r25</span><span class='hljl-t'>
 </span><span class='hljl-n'>dy</span><span class='hljl-p'>[</span><span class='hljl-ni'>20</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>r25</span><span class='hljl-oB'>+</span><span class='hljl-n'>r24</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span>
</pre>


<pre class="output">
f &#40;generic function with 1 method&#41;
</pre>



<pre class='hljl'>
<span class='hljl-n'>u0</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>zeros</span><span class='hljl-p'>(</span><span class='hljl-ni'>20</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>u0</span><span class='hljl-p'>[</span><span class='hljl-ni'>2</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.2</span><span class='hljl-t'>
</span><span class='hljl-n'>u0</span><span class='hljl-p'>[</span><span class='hljl-ni'>4</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.04</span><span class='hljl-t'>
</span><span class='hljl-n'>u0</span><span class='hljl-p'>[</span><span class='hljl-ni'>7</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.1</span><span class='hljl-t'>
</span><span class='hljl-n'>u0</span><span class='hljl-p'>[</span><span class='hljl-ni'>8</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.3</span><span class='hljl-t'>
</span><span class='hljl-n'>u0</span><span class='hljl-p'>[</span><span class='hljl-ni'>9</span><span class='hljl-p'>]</span><span class='hljl-t'>  </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.01</span><span class='hljl-t'>
</span><span class='hljl-n'>u0</span><span class='hljl-p'>[</span><span class='hljl-ni'>17</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.007</span><span class='hljl-t'>
</span><span class='hljl-n'>prob</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ODEProblem</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-n'>u0</span><span class='hljl-p'>,(</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-nfB'>60.0</span><span class='hljl-p'>),</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sol</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>solve</span><span class='hljl-p'>(</span><span class='hljl-n'>prob</span><span class='hljl-p'>,</span><span class='hljl-nf'>Rodas5</span><span class='hljl-p'>())</span>
</pre>


<pre class="output">
retcode: Success
Interpolation: 3rd order Hermite
t: 29-element Array&#123;Float64,1&#125;:
  0.0
  0.0013845590497824308
  0.003242540880561935
  0.007901604453787887
  0.016011571579446486
  0.02740615501991768
  0.04461209330420329
  0.07720628505724063
  0.11607397409282896
  0.1763063424034578
  ⋮
  3.67632101025685
  5.344944969239927
  7.985768407909472
 11.981250465617169
 17.4525037309389
 24.88224633682058
 34.664620859612675
 47.27762914412602
 60.0
u: 29-element Array&#123;Array&#123;Float64,1&#125;,1&#125;:
 &#91;0.0, 0.2, 0.0, 0.04, 0.0, 0.0, 0.1, 0.3, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0, 0
.0, 0.0, 0.007, 0.0, 0.0, 0.0&#93;
 &#91;0.0002935083676916062, 0.19970649101355778, 1.693835912963263e-10, 0.0397
0673366392418, 1.1424841091201038e-7, 1.0937647553427759e-7, 0.099999666764
40009, 0.3000003350396963, 0.00999998209858389, 5.6037724774041996e-9, 6.04
7938337401518e-9, 1.004757080331196e-8, 5.9812281205517996e-12, 6.239553394
674922e-9, 2.3022538431539757e-10, 3.129330582279999e-17, 0.006999999417662
247, 5.823377531818463e-10, 3.8240538930006997e-10, 6.930819262664895e-14&#93;
 &#91;0.0006836584262738185, 0.19931633629695783, 1.9606016564783855e-10, 0.039
31745222652303, 2.0821314359971526e-7, 2.5242637698192594e-7, 0.09999884087
661927, 0.30000116344748634, 0.009999897466494299, 2.060514512836631e-8, 1.
6844788957556812e-8, 8.136618249763386e-8, 1.0724538885842355e-10, 6.486750
944005252e-8, 3.1010416395449586e-9, 3.098650817258985e-17, 0.0069999964441
40231, 3.555859768109285e-9, 2.064386087388623e-9, 2.0476275988812504e-12&#93;
 &#91;0.0016428178251749216, 0.1983571313918958, 2.6245923327739e-10, 0.0383623
1113129213, 3.515272953820738e-7, 4.664572294956194e-7, 0.09999544921432231
, 0.30000456320447505, 0.009999473651103152, 4.4964381405774644e-8, 3.33502
8282086443e-8, 4.812857607411457e-7, 1.44096025408721e-9, 4.444463679809393
e-7, 3.7457507792328246e-8, 3.023375117340986e-17, 0.006999981334745823, 1.
8665254176770853e-8, 1.1746740499997834e-8, 6.886036636810248e-11&#93;
 &#91;0.0032462847242510895, 0.19675344603716596, 3.7349153139939496e-10, 0.036
76859844153611, 4.3203353387790993e-7, 5.784968999241242e-7, 0.099987536470
72418, 0.30001250073127955, 0.0099984127906682, 5.847792637527643e-8, 4.228
372979822697e-8, 1.5155868686349532e-6, 8.524978790566793e-9, 1.46153444758
63662e-6, 2.1422142790901277e-7, 2.8977729013034015e-17, 0.0069999433444137
9, 5.66555862100937e-8, 4.4368490665966554e-8, 1.0618428053776192e-9&#93;
 &#91;0.005361392688498683, 0.19463777435511, 5.199898634456225e-10, 0.03466792
339483551, 4.3013501413959053e-7, 5.629734362718953e-7, 0.09997580963070891
, 0.30002429018620363, 0.009996820355063572, 5.777644600442459e-8, 4.151145
8414889325e-8, 3.075243081277304e-6, 2.7267220452845075e-8, 2.9888960752504
05e-6, 6.762048663689083e-7, 2.7322164351120048e-17, 0.006999886274028278, 
1.1372597172141333e-7, 1.1359724052467222e-7, 7.943531994069594e-9&#93;
 &#91;0.008274965064030285, 0.19172294835291767, 7.218269648981172e-10, 0.03177
417980097769, 3.98260843239196e-7, 5.005231442113355e-7, 0.0999593539433593
2, 0.3000408990194988, 0.009994606735402435, 5.178145361951533e-8, 3.714992
467251187e-8, 5.229419114367511e-6, 6.871429897641355e-8, 5.040637321652717
5e-6, 1.6909615266613677e-6, 2.5041573813257427e-17, 0.006999806991156842, 
1.9300884315678748e-7, 2.3808904098709325e-7, 4.440909272869548e-8&#93;
 &#91;0.01300266954290988, 0.1869919466035047, 1.0493234388726091e-9, 0.0270783
35801454113, 3.596574054568645e-7, 4.230621124193838e-7, 0.0999319172001279
2, 0.30006878678581594, 0.009990985413337174, 4.4266325041299826e-8, 3.1717
471243607364e-8, 8.707147429503103e-6, 1.7539879080367005e-7, 8.15954110305
292e-6, 4.282354221742325e-6, 2.134072850846997e-17, 0.006999677462815518, 
3.2253718448129225e-7, 4.2925309447948683e-7, 2.484237391964801e-7&#93;
 &#91;0.017549505096292266, 0.18244011104288263, 1.3641805131831696e-9, 0.02256
3571980700702, 3.361162857847711e-7, 3.7150094810967537e-7, 0.0999030946371
2853, 0.300098365438841, 0.009987255300925597, 3.9316875533248874e-8, 2.813
081375953118e-8, 1.2231148883147214e-5, 3.34623584761564e-7, 1.103334762344
7913e-5, 8.113768214597823e-6, 1.778259444568051e-17, 0.006999544244978277,
 4.557550217235584e-7, 5.171121300549293e-7, 7.09178447867133e-7&#93;
 &#91;0.0228606589237139, 0.17711993568649875, 1.731845524723872e-9, 0.01729532
122114822, 3.1879076031327097e-7, 3.280176124825831e-7, 0.09986309179339616
, 0.300139895927231, 0.009982163707599944, 3.519652846606691e-8, 2.51367691
36586453e-8, 1.6956870484251417e-5, 6.253074103684872e-7, 1.439190880113898
3e-5, 1.503681700195539e-5, 1.3630629199414105e-17, 0.006999362663103371, 6
.373368966290782e-7, 5.039634221292022e-7, 1.619650976481737e-6&#93;
 ⋮
 &#91;0.03892999698723858, 0.16044023738492405, 2.851149270321245e-9, 0.0031795
600490061044, 3.0972898091824323e-7, 2.602522850073365e-7, 0.09808114628680
965, 0.30211293726179844, 0.00975642790788472, 2.885705502284382e-8, 2.0507
65092305705e-8, 0.00021451968557192194, 2.4313107185131944e-5, 2.9844863323
202143e-5, 0.0005806905788187015, 2.5058455689081163e-18, 0.006991260060189
278, 8.739939810724537e-6, 5.644795247306144e-7, 1.2098731154429619e-5&#93;
 &#91;0.039731796047813156, 0.1593487586247586, 2.910031089549733e-9, 0.0032678
396978631824, 3.0340454783426414e-7, 2.533682739802742e-7, 0.09727159545005
896, 0.30301666794722515, 0.009654512439051292, 2.803900061546711e-8, 1.991
4463772627238e-8, 0.000303486670316648, 3.515693085822018e-5, 2.88433246518
7963e-5, 0.0008554036328153492, 2.5754197123442302e-18, 0.00698754630983515
65, 1.245369016484512e-5, 6.384025408091165e-7, 1.4123180606949498e-5&#93;
 &#91;0.04094318642892073, 0.15768509237182085, 2.998971461318937e-9, 0.0034038
15741811088, 2.9398821558655013e-7, 2.4320600465068163e-7, 0.09603188750370
362, 0.3043989905948464, 0.009500566586930251, 2.684636982205443e-8, 1.9049
770161614344e-8, 0.0004379930491040918, 5.130947934360084e-5, 2.73117738878
48295e-5, 0.0012863758134774614, 2.6825838991980204e-18, 0.0069818694792832
51, 1.8130520716750667e-5, 7.249055743597863e-7, 1.6655500431521582e-5&#93;
 &#91;0.042654560651108626, 0.1553017077035719, 3.124615480512117e-9, 0.0036015
751816334397, 2.8144780078033687e-7, 2.2972055207124975e-7, 0.0942434704647
6046, 0.30638970261426646, 0.009282897888913216, 2.528730372862258e-8, 1.79
19685179564033e-8, 0.000628530891181199, 7.356402437553935e-5, 2.5311076540
045087e-5, 0.0019296589138284859, 2.8384402467274616e-18, 0.006973700751202
073, 2.6299248797929046e-5, 8.253072035773858e-7, 1.9841699955945823e-5&#93;
 &#91;0.04479393132409089, 0.1522628962759311, 3.281711110095147e-9, 0.00385892
21418967377, 2.6689444184296353e-7, 2.1414676097184994e-7, 0.09194251773076
727, 0.3089455003617426, 0.009010154457710506, 2.3520931505423736e-8, 1.663
978888926712e-8, 0.0008681120475472489, 0.00010022358375904649, 2.305697649
42253e-5, 0.0027942301640086443, 3.0412581618190092e-18, 0.0069632199666104
46, 3.6780033389556585e-5, 9.439336027242336e-7, 2.3887359303821783e-5&#93;
 &#91;0.04738795680443153, 0.14848196867160532, 3.4722646313007398e-9, 0.004187
711078360779, 2.506414878535045e-7, 1.9688141906496757e-7, 0.08904810579053
157, 0.3121528524878342, 0.008678047289621835, 2.1605580834173356e-8, 1.525
2579587683115e-8, 0.0011615457673835054, 0.00013035939033899277, 2.06304835
8889216e-5, 0.003939836161922413, 3.3003802689174486e-18, 0.006950068617092
346, 4.993138290765649e-5, 1.096031158215806e-6, 2.939147027175925e-5&#93;
 &#91;0.05036882011197725, 0.14399340521526238, 3.691354467227887e-9, 0.0045913
22401963945, 2.3347510372053892e-7, 1.788424393930541e-7, 0.085567494759836
05, 0.31600005128850905, 0.008293774660545702, 1.964965282178699e-8, 1.3836
778039561243e-8, 0.0015041058356833838, 0.0001612874061394357, 1.8176659392
1036e-5, 0.0054015179236651485, 3.618470706343953e-18, 0.006934279702559454
, 6.572029744054831e-5, 1.2906023005324613e-6, 3.683937032755865e-5&#93;
 &#91;0.053649127243805966, 0.13885256756140255, 3.932627684418802e-9, 0.005072
852275710312, 2.1601059359540682e-7, 1.6077417122164273e-7, 0.0815216578847
2086, 0.3204605691433712, 0.007866252082050794, 1.772858492881216e-8, 1.244
7110598095664e-8, 0.0018899285762960995, 0.00018974902469467097, 1.57995593
14561062e-5, 0.007213605381582541, 3.997969593556813e-18, 0.006915930520936
276, 8.406947906372734e-5, 1.5343173253873047e-6, 4.6708235594326756e-5&#93;
 &#91;0.05646254720111089, 0.1342484201017065, 4.1397337819370874e-9, 0.0055231
39725449666, 2.0189808613242023e-7, 1.4645450805872716e-7, 0.07784249512128
202, 0.32450753046673436, 0.007494014171885585, 1.6222966604690566e-8, 1.13
58663790555432e-8, 0.002230505164655582, 0.00020871630726228214, 1.39693487
37895171e-5, 0.008964884997558377, 4.3528459894333626e-18, 0.00689921972266
9937, 0.00010078027733006581, 1.7721521050391386e-6, 5.682962012835069e-5&#93;
</pre>



<pre class='hljl'>
<span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>)</span>
</pre>


<img src=""  />


<pre class='hljl'>
<span class='hljl-nf'>plot</span><span class='hljl-p'>(</span><span class='hljl-n'>sol</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>xscale</span><span class='hljl-oB'>=:</span><span class='hljl-n'>log10</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>tspan</span><span class='hljl-oB'>=</span><span class='hljl-p'>(</span><span class='hljl-nfB'>1e-6</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-ni'>60</span><span class='hljl-p'>),</span><span class='hljl-t'> </span><span class='hljl-n'>layout</span><span class='hljl-oB'>=</span><span class='hljl-p'>(</span><span class='hljl-ni'>3</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>))</span>
</pre>


<img src=""  />

<h2>Geometric Properties</h2>
<h3>Linear Ordinary Differential Equations</h3>
<p>The simplest ordinary differential equation is the scalar linear ODE, which is given in the form</p>
<p class="math">\[
u' = \alpha u
\]</p>
<p>We can solve this by noticing that <span class="math">$(e^{\alpha t})^\prime = \alpha e^{\alpha t}$</span> satisfies the differential equation and thus the general solution is:</p>
<p class="math">\[
u(t) = u(0)e^{\alpha t}
\]</p>
<p>From the analytical solution we have that:</p>
<ul>
<li><p>If <span class="math">$Re(\alpha) > 0$</span> then <span class="math">$u(t) \rightarrow \infty$</span> as <span class="math">$t \rightarrow \infty$</span></p>
</li>
<li><p>If <span class="math">$Re(\alpha) < 0$</span> then <span class="math">$u(t) \rightarrow 0$</span> as <span class="math">$t \rightarrow \infty$</span></p>
</li>
<li><p>If <span class="math">$Re(\alpha) = 0$</span> then <span class="math">$u(t)$</span> has a constant or periodic solution.</p>
</li>
</ul>
<p>This theory can then be extended to multivariable systems in the same way as the discrete dynamics case. Let <span class="math">$u$</span> be a vector and have</p>
<p class="math">\[
u' = Au
\]</p>
<p>be a linear ordinary differential equation. Assuming <span class="math">$A$</span> is diagonaliziable, we diagonalize <span class="math">$A = P^{-1}DP$</span> to get</p>
<p class="math">\[
Pu' = DPu
\]</p>
<p>and change coordinates <span class="math">$z = Pu$</span> so that we have</p>
<p class="math">\[
z' = Dz
\]</p>
<p>which decouples the equation into a system of linear ordinary differential equations which we solve individually. Thus we see that, similarly to the discrete dynamical system, we have that:</p>
<ul>
<li><p>If all of the eigenvalues negative, then <span class="math">$u(t) \rightarrow 0$</span> as <span class="math">$t \rightarrow \infty$</span></p>
</li>
<li><p>If any eigenvalue is positive, then <span class="math">$u(t) \rightarrow \infty$</span> as <span class="math">$t \rightarrow \infty$</span></p>
</li>
</ul>
<h3>Nonlinear Ordinary Differential Equations</h3>
<p>As with discerte dynamical systems, the geometric properties extend locally to the linearization of the continuous dynamical system as defined by:</p>
<p class="math">\[
u' = \frac{df}{du} u
\]</p>
<p>where <span class="math">$\frac{df}{du}$</span> is the Jacobian of the system. This is a consequence of the Hartman-Grubman Theorem.</p>
<h2>Numerically Solving Ordinary Differential Equations</h2>
<h3>Euler&#39;s Method</h3>
<p>To numerically solve an ordinary differential equation, one turns the continuous equation into a discrete equation by <em>discretizing</em> it. The simplest discretization is the <em>Euler method</em>. The Euler method can be thought of as a simple approximation replacing <span class="math">$dt$</span> with a small non-infinitesimal <span class="math">$\Delta t$</span>. Thus we can approximate</p>
<p class="math">\[
f(u,p,t) = u' = \frac{du}{dt} \approx \frac{\Delta u}{\Delta t}
\]</p>
<p>and now since <span class="math">$\Delta u = u_{n+1} - u_n$</span> we have that</p>
<p class="math">\[
\Delta t f(u,p,t) = u_{n+1} - u_n
\]</p>
<p>We need to make a choice as to where we evaluate <span class="math">$f$</span> at. The simplest approximation is to evaluate it at <span class="math">$t_n$</span> with <span class="math">$u_n$</span> where we already have the data, and thus we re-arrange to get</p>
<p class="math">\[
u_{n+1} = u_n + \Delta t f(u,p,t)
\]</p>
<p>This is the Euler method.</p>
<p>We can interpret it more rigorously by looking at the Taylor series expansion. First write out the Taylor series for the ODE&#39;s solution in the near future:</p>
<p class="math">\[
u(t+\Delta t) = u(t) + \Delta t u'(t) + \frac{\Delta t^2}{2} u''(t) + \ldots
\]</p>
<p>Recall that <span class="math">$u' = f(u,p,t)$</span> by the definition of the ODE system, and thus we have that</p>
<p class="math">\[
u(t+\Delta t) = u(t) + \Delta t f(u,p,t) + \mathcal{O}(\Delta t^2)
\]</p>
<p>This is a first order approximation because the error in our step can be expresed as an error in the derivative, i.e.</p>
<p class="math">\[
\frac{u(t + \Delta t) - u(t)}{\Delta t} = f(u,p,t) + \mathcal{O}(\Delta t)
\]</p>
<h3>Higher Order Methods</h3>
<p>We can use this analysis to extend our methods to higher order approximation by simply matching the Taylor series to a higher order. Intuitively, when we developed the Euler method we had to make a choice:</p>
<p class="math">\[
u_{n+1} = u_n + \Delta t f(u,p,t)
\]</p>
<p>where do we evaluate <span class="math">$f$</span>? One may think that the best derivative approximation my come from the middle of the interval, in which case we might want to evaluate it at <span class="math">$t + \frac{\Delta t}{2}$</span>. To do so, we can use the Euler method to approximate the value at <span class="math">$t + \frac{\Delta t}{2}$</span> and then use that value to approximate the derivative at <span class="math">$t + \frac{\Delta t}{2}$</span>. This looks like:</p>
<p class="math">\[
k_1 = f(u_n,p,t)\\
k_2 = f(u_n + \frac{\Delta t}{2} k_1,p,t + \frac{\Delta t}{2})\\
u_{n+1} = u_n + \Delta t k_2
\]</p>
<p>which we can also write as:</p>
<p class="math">\[
u_{n+1} = u_n + \Delta t f(t + \frac{\Delta t}{2},u_n + \frac{\Delta t}{2} f_n)
\]</p>
<p>where <span class="math">$f_n = f(u_n,p,t)$</span>. If we do the two-dimensional Taylor expansion we get:</p>
<p class="math">\[
u_{n+1} = u_n + \Delta t f_n + \frac{\Delta t^2}{2}(f_t + f_u f)(u_n,p,t)\\
+ \frac{\Delta t^3}{6} (f_{tt} + 2f_{tu}f + f_{uu}f^2)(u_n,p,t)
\]</p>
<p>which when we compare against the true Taylor series:</p>
<p class="math">\[
u(t+\Delta t) = u_n + \Delta t f(u_n,p,t) + \frac{\Delta t^2}{2}(f_t + f_u f)(u_n,p,t)
+ \frac{\Delta t^3}{6}(f_{tt} + 2f_{tu} + f_{uu}f^2 + f_t f_u + f_u^2 f)(u_n,p,t)
\]</p>
<p>and thus we see that</p>
<p class="math">\[
u(t + \Delta t) - u_n = \mathcal{O}(\Delta t^3)
\]</p>
<h3>Runge-Kutta Methods</h3>
<p>More generally, Runge-Kutta methods are of the form:</p>
<p class="math">\[
k_1 = f(u_n,p,t)\\
k_2 = f(u_n + \Delta t (a_{21} k_1),p,t + \Delta t c_1)\\
k_3 = f(u_n + \Delta t (a_{31} k_1 + a_{32} k_2),p,t + \Delta t c_2)\\
\vdots \\
u_{n+1} = u_n + \Delta t (b_1 k_1 + \ldots + b_s k_s)
\]</p>
<p>where <span class="math">$s$</span> is the number of stages. These tableaus can be expressed as a tableau:</p>
<p><img src="https://en.wikipedia.org/wiki/List_of_Runge&#37;E2&#37;80&#37;93Kutta_methods" alt="" /></p>
<p>The order of the Runge-Kutta method is simply the number of terms in the Taylor series that ends up being matched by the resulting expansion. For example, for the 4th order you can expand out and see that the following equations need to be satisfied:</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95117136-105ae780-0716-11eb-9f6a-49fecf7adbeb.PNG" alt="" /></p>
<p>The classic Runge-Kutta method is also known as RK4 and is the following 4th order method:</p>
<p class="math">\[
k_1 = f(u_n,p,t)\\
k_2 = f(u_n + \frac{\Delta t}{2} k_1,p,t + \frac{\Delta t}{2})\\
k_3 = f(u_n + \frac{\Delta t}{2} k_2,p,t + \frac{\Delta t}{2})\\
k_4 = f(u_n + \Delta t k_3,p,t + \Delta t)\\
u_{n+1} = u_n + \frac{1}{6}(k_1 + 2 k_2 + 2 k_3 + k_4)\\
\]</p>
<p>While it&#39;s widely known and simple to remember, it&#39;s not necessarily good. The way to judge a Runge-Kutta method is by looking at the size of the coefficient of the next term in the Taylor series: if it&#39;s large then the true error can be larger, even if it matches another one asymtopically.</p>
<h2>What Makes a Good Method?</h2>
<h3>Leading Truncation Coeffcients</h3>
<p>For given orders of explicit Runge-Kutta methods, lower bounds for the number of <code>f</code> evaluations &#40;stages&#41; required to receive a given order are known:</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95117078-f8836380-0715-11eb-9acf-0626338307d1.PNG" alt="" /></p>
<p>While uninuitive, using the method is not necessarily the one that reduces the coefficient the most. The reason is because what is attempted in ODE solving is precisely the opposite of the analysis. In the ODE analysis, we&#39;re looking at behavior as <span class="math">$\Delta t \rightarrow 0$</span>. However, when efficiently solving ODEs, we want to use the largest <span class="math">$\Delta t$</span> which satisfies error tolerances.</p>
<p>The most widely used method is the Dormand-Prince 5th order Runge-Kutta method, whose tableau is represented as:</p>
<p><img src="http://rotordynamics.files.wordpress.com/2014/05/new-picture6.png" alt="" /></p>
<p>Notice that this method takes 7 calls to <code>f</code> for 5th order. The key to this method is that it has optimized leading truncation error coefficients, under some extra assumptions which allow for the analysis to be simplified.</p>
<h3>Looking at the Effects of RK Method Choices and Code Optimizations</h3>
<p>Pulling from the <a href="https://github.com/SciML/SciMLBenchmarks.jl">SciML Benchmarks</a>, we can see the general effect of these different properties on a given set of Runge-Kutta methods:</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95118000-7c8a1b00-0717-11eb-8080-2179da500cd2.PNG" alt="" /></p>
<p>Here, the order of the method is given in the name. We can see one immediate factor is that, as the requested error in the calculation decreases, the higher order methods become more efficient. This is because to decrease error, you decrease <span class="math">$\Delta t$</span>, and thus the exponent difference with respect to <span class="math">$\Delta t$</span> has more of a chance to pay off for the extra calls to <code>f</code>. Additionally, we can see that order is not the only determining factor for efficiency: the Vern8 method seems to have a clear approximate 2.5x performance advantage over the whole span of the benchmark compared to the DP8 method, even though both are 8th order methods. This is because of the leading truncation terms: with a small enough <span class="math">$\Delta t$</span>, the more optimized method &#40;Vern8&#41; will generally have low error in a step for the same <span class="math">$\Delta t$</span> because the coefficients in the expansion are generally smaller.</p>
<p>This is a factor which is generally ignored in high level discussions of numerical differential equations, but can lead to orders of magnitude differences&#33; This is highlighted in the following plot:</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95118457-544eec00-0718-11eb-8c19-f402e2cb8842.PNG" alt="" /></p>
<p>Here we see ODEInterface.jl&#39;s ODEInterfaceDiffEq.jl wrapper into the SciML common interface for the standard <code>dopri</code> method from Fortran, and ODE.jl, the original ODE solvers in Julia, have a performance disadvantage compared to the DifferentialEquations.jl methods due in part to some of the coding performance pieces that we discussed in the first few lectures.</p>
<p>Specifically, a large part of this can be attributed to inlining of the higher order functions, i.e. ODEs are defined by a user function and then have to be called from the solver. If the solver code is compiled as a shared library ahead of time, like is commonly done in C&#43;&#43; or Fortran, then there can be a function call overhead that is eliminted by JIT compilation optimizing across the function call barriers &#40;known as interprocedural optimization&#41;. This is one way which a JIT system can outperform an AOT &#40;ahead of time&#41; compiled system in real-world code &#40;for completeness, two other ways are by doing full function specialization, which is something that is <a href="https://scalac.io/specialized-generics-object-instantiation/">not generally possible in AOT languages given that you cannot know all types ahead of time for a fully generic function</a>, and <a href="https://github.com/dyu/ffi-overhead#results-500m-calls">calling C itself, i.e. c-ffi &#40;foreign function interface&#41;, can be optimized using the runtime information of the JIT compiler to outperform C&#33;</a>&#41;.</p>
<p>The other performance difference being shown here is due to optimization of the method. While a slightly different order, we can see a clear difference in the performance of RK4 vs the coefficient optimized methods. It&#39;s about the same order of magnitude as &quot;highly optimized code differences&quot;, showing that what&#39;s different about a both the Runge-Kutta coefficients and the code implementation can both have a significant impact on performance.</p>
<p>Taking a look at what happens when interpreted languages get involved highlights some of the code challenges in this domain. Let&#39;s take a look at for example the results when simulating 3 ODE systems with the various RK methods:</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95131785-b1549d00-072c-11eb-8d2a-490f69a4b99f.PNG" alt="" /></p>
<p>We see that using interpreted languages introduces around a 50x-100x performance penalty. If the you recall in your previous lecture, the discrete dynamical system that was being simulated was the 3-dimensional Lorenz equation discretized by Euler&#39;s method, meaning that the performance of that implementation is a good proxy for understanding the performance differences in this graph. Recall that in previous lectures we saw an approximately 5x performance advantage when specializing on the system function and size and around 10x by reducing allocations: these features account for the performance differences noticed between library implementations, which are then compounded by the use of different RK methods &#40;note that R uses &quot;call by copy&quot; which even further increases the memory usages and makes standard usage of the language incompatible with mutating function calls&#33;&#41;.</p>
<h3>Stability of a Method</h3>
<p>Simply having an order on the truncation error does not imply convergence of the method. The disconnect is that the errors at a given time point may not dissipate. What also needs to be checked is the asymtopic behavior of a disturbance. To see this, one can utilize the linear test problem:</p>
<p class="math">\[
u' = \alpha u
\]</p>
<p>and ask the question, does the discrete dynamical system defined by the discretized ODE end up going to zero? You would hope that the discretized dynamical system and the continuous dynamical system have the same properties in this simple case, and this is known as linear stability analysis of the method.</p>
<p>As an example, take a look at the Euler method. Recall that the Euler method was given by:</p>
<p class="math">\[
u_{n+1} = u_n + \Delta t f(u_n,p,t)
\]</p>
<p>When we plug in the linear test equation, we get that</p>
<p class="math">\[
u_{n+1} = u_n + \Delta t \alpha u_n
\]</p>
<p>If we let <span class="math">$z = \Delta t \alpha$</span>, then we get the following:</p>
<p class="math">\[
u_{n+1} = u_n + z u_n = (1+z)u_n
\]</p>
<p>which is stable when <span class="math">$z$</span> is in the shifted unit circle. This means that, as a necessary condition, the step size <span class="math">$\Delta t$</span> needs to be small enough that <span class="math">$z$</span> satisfies this condition, placing a stepsize limit on the method.</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95117231-3c766880-0716-11eb-9069-039253bcebda.PNG" alt="" /></p>
<p>If <span class="math">$\Delta t$</span> is ever too large, it will cause the equation to overshoot zero, which then causes oscillations that spiral out to infinity.</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95132604-0d6bf100-072e-11eb-8af5-663512a0db14.PNG" alt="" /></p>
<p><img src="https://user-images.githubusercontent.com/1814174/95132963-9125dd80-072e-11eb-878e-61f77a20d03e.gif" alt="" /></p>
<p>Thus the stability condition places a hard constraint on the allowed <span class="math">$\Delta t$</span> which will result in a realistic simulation.</p>
<p>For reference, the stability regions of the 2nd and 4th order Runge-Kutta methods that we discussed are as follows:</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95117286-56b04680-0716-11eb-9c6a-07fc4d190a09.PNG" alt="" /></p>
<h3>Interpretation of the Linear Stability Condition</h3>
<p>To interpret the linear stability condition, recall that the linearization of a system interprets the dynamics as locally being due to the Jacobian of the system. Thus</p>
<p class="math">\[
u' = f(u,p,t)
\]</p>
<p>is locally equivalent to</p>
<p class="math">\[
u' = \frac{df}{du}u
\]</p>
<p>You can understand the local behavior through diagonalizing this matrix. Therefore, the scalar for the linear stability analysis is performing an analysis on the eigenvalues of the Jacobian. The method will be stable if the largest eigenvalues of df/du are all within the stability limit. This means that stability effects are different throughout the solution of a nonlinear equation and are generally understood locally &#40;though different more comprehensive stability conditions exist&#33;&#41;.</p>
<h3>Implicit Methods</h3>
<p>If instead of the Euler method we defined <span class="math">$f$</span> to be evaluated at the future point, we would receive a method like:</p>
<p class="math">\[
u_{n+1} = u_n + \Delta t f(u_{n+1},p,t+\Delta t)
\]</p>
<p>in which case, for the stability calculation we would have that</p>
<p class="math">\[
u_{n+1} = u_n + \Delta t \alpha u_n
\]</p>
<p>or</p>
<p class="math">\[
(1-z) u_{n+1} = u_n
\]</p>
<p>which means that</p>
<p class="math">\[
u_{n+1} = \frac{1}{1-z} u_n
\]</p>
<p>which is stable for all <span class="math">$Re(z) < 0$</span> a property which is known as A-stability. It is also stable as <span class="math">$z \rightarrow \infty$</span>, a property known as L-stability. This means that for equations with very ill-conditioned Jacobians, this method is still able to be use reasonably large stepsizes and can thus be efficient.</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95117191-28326b80-0716-11eb-8e17-889308bdff53.PNG" alt="" /></p>
<h3>Stiffness and Timescale Separation</h3>
<p>From this we see that there is a maximal stepsize whenever the eigenvalues of the Jacobian are sufficiently large. It turns out that&#39;s not an issue if the phonomena we fast to see is fast, since then the total integration time tends to be small. However, is we have some equations with both fast modes and slow modes, like the Robertson equation, then it is very difficult because in order to resolve the slow dynamics over a long timespan, one needs to ensure that the fast dynamics do not diverge. This is a property known as stiffness. Stiffness can thus be approximated in some sense by the condition number of the Jacobian. The condition number of a matrix is its maximal eigenvalue divided by its minimal eigenvalue and gives an rough measure of the local timescale separations. If this value is large and one wants to resolve the slow dynamics, then explict integrators, like the explicit Runge-Kutta methods described before, have issues with stability. In this case implicit integrators &#40;or other forms of stabilized stepping&#41; are required in order to efficiently reach the end time step.</p>
<p><img src="https://user-images.githubusercontent.com/1814174/95132552-f6c59a00-072d-11eb-881e-24364b7b728f.PNG" alt="" /></p>
<h2>Exploiting Continuity</h2>
<p>So far, we have looked at ordinary differential equations as a <span class="math">$\Delta t \rightarrow 0$</span> formulation of a discrete dynamical system. However, continuous dynamics and discrete dynamics have very different characteristics which can be utilized in order to arrive at simpler models and faster computations.</p>
<h3>Geometric Properties: No Jumping and the Poincaré–Bendixson theorem</h3>
<p>In terms of geometric properties, continuity places a large constraint on the possible dynamics. This is because of the physical constraint on &quot;jumping&quot;, i.e. flows of differential equations cannot jump over each other. If you are ever at some point in phase space and <span class="math">$f$</span> is not explicitly time-dependent, then the direction of <span class="math">$u'$</span> is uniquely determined &#40;given reasonable assumptions on <span class="math">$f$</span>&#41;, meaning that flow lines &#40;solutions to the differential equation&#41; can never cross.</p>
<p>A result from this is the Poincaré–Bendixson theorem, which states that, with any arbitrary &#40;but nice&#41; two dimensional continuous system, you can only have 3 behaviors:</p>
<ul>
<li><p>Steady state behavior</p>
</li>
<li><p>Divergence</p>
</li>
<li><p>Periodic orbits</p>
</li>
</ul>
<p>A simple proof by picture shows this.</p>


        <HR/>
        <div class="footer">
          <p>
            Published from <a href="discretizing_odes.jmd">discretizing_odes.jmd</a>
            using <a href="http://github.com/JunoLab/Weave.jl">Weave.jl</a> v0.10.3 on 2020-10-08.
          </p>
        </div>
      </div>
    </div>
  </div>
</BODY>

</HTML>
